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Abstract 

Differential quantities, including normals, curvatures, principal di- 
rections, and associated matrices, play a fundamental role in geo- 
metric processing and physics-based modeling. Computing these 
differential quantities consistently on surface meshes is important 
and challenging, and some existing methods often produce incon- 
sistent results and require ad hoc fixes. In this paper, we show that 
the computation of the gradient and Hessian of a height function 
provides the foundation for consistently computing the differential 
quantities. We derive simple, explicit formulas for the transforma- 
tions between the first- and second-order differential quantities (i.e., 
normal vector and principal curvature tensor) of a smooth surface 
and the first- and second-order derivatives (i.e., gradient and Hes- 
sian) of its corresponding height function. We then investigate a 
general, flexible numerical framework to estimate the derivatives of 
the height function based on local polynomial fittings formulated 
as weighted least squares approximations. We also propose an it- 
erative fitting scheme to improve accuracy. This framework gen- 
eralizes polynomial fitting and addresses some of its accuracy and 
stability issues, as demonstrated by our theoretical analysis as well 
as experimental results. 
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ometry and Object Modeling — Algorithms, Design, Experimenta- 
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1 Introduction 

Computing normals and curvatures is a fundamental problem for 
many geometric and numerical computations, including feature de- 
tection, shape retrieval, shape registration or matching, surface 
fairing, surface mesh adaptation or remeshing, front tracking and 
moving meshes. In recent years, a number of methods have 
been introduced for the computation of the differential quantities 
(see e.g., pj iubin 1995l lMeek and Walton 2000l|Meyer et al. 2002| 
|Cazals ancTP ouget 2005 1). However, some of the existing methods 
may produce inconsistent results. For example, when estimating 
the mean curvature using the cotangent formula and estimating the 
Gaussian curvature using the angle deficit | Meyer et al. 2002), the 
principal curvatures obtained from these mean and Gaussian curva- 
tures are not guaranteed to be real numbers. Such inconsistencies 
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often require ad hoc fixes to avoid crashing of the code, and their 
effects on the accuracy of the applications are difficult to analyze. 

The ultimate goal of this work is to investigate a mathematically 
sound framework that can compute the differential quantities con- 
sistently (i.e., satisfying the intrinsic constraints) with provable 
convergence on general surface meshes, while being flexible and 
easy to implement. This is undoubtly an ambitious goal. Al- 
though we may have not fully achieved the goal, we make some 
contributions toward it. First, using the singular value decompo- 
sition |Golub and Van Loan 1996 1 of the Jacobian matrix, we de- 
rive explicit formulas for the transformations between the first- and 
second-order differential quantities of a smooth surface (i.e., nor- 
mal vector and principal curvature tensor) and the first- and second- 
order derivatives of its corresponding height function (i.e., gradient 
and Hessian). We also give the explicit formulas for the transforma- 
tions of the gradient and Hessian under a rotation of the coordinate 
system. These transformations can be obtained without forming the 
shape operator and the associated computation of its eigenvalues or 
eigenvectors. We then reduce the problem, in a consistent man- 
ner, to the computation of the gradient and Hessian of a function 
in two dimensions, which can then be computed more readily by 
customizing classical numerical techniques. 

Our second contribution is to develop a general and flexible 
method for differentiating a height function, which can be viewed 
as a generalization of the polynomial fitting of osculating jets 
(Cazals and Pouget 2005 ]. Our method is based on the Taylor se- 
ries expansions of a function and its derivatives, then solved by a 
weighted least squares formulation. We also propose an iterative- 
fitting method that improves the accuracy of fittings. Further- 
more, we address the numerical stability issues by a systematic 
point-selection strategy and a numerical solver with safeguard. We 
present experimental results to demonstrate the accuracy and stabil- 
ity of our approach for fittings up to degree six. 

The remainder of this paper is organized as follows. Section|2]re- 
views some related work on the computation of differential quanti- 
ties of discrete surfaces. Section[5]analyzes the stability and consis- 
tency of classical formulas for computing differential quantities for 
continuous surfaces and establishes the theoretical foundation for 
consistent computations using a height function. Section|4]presents 
a general framework and a unified analysis for computing the gra- 
dient and Hessian of a height function based on polynomial fittings, 
including an iterative-fitting scheme. Section |3]presents some ex- 
perimental results to demonstrate the accuracy and stability of our 
approach. Section|6]concludes the paper with a discussion on future 
research directions. 

2 Related Work 

Many methods have been proposed for computing or estimating the 
first- and second-order differential quantities of a surface. In recent 
years, there have been significant interests in the convergence and 
consistency of these methods. We do not attempt to give a compre- 
hensive review of these methods but consider only a few of them 
that are more relevant to our proposed approach; readers are re- 
ferred to |Petitjean 2002) and IGatzke and G rimm 2006] for com- 



prehensive surveys. Among the existing methods, many of them 
estimate the different quantities separately of each other. For the 
estimation of the normals, a common practice is to estimate ver- 
tex normals using a weighted average of face normals, such as area 
weighting or angle weighting. These methods in general are only 
first-order accurate, although they are the most efficient. 



For curvature estimation, Meyer et al. 1 2002 1 proposed some 
discrete operators for computing normals and curvatures, among 
which the more popular one is the cotangent formula for mean cur- 
vatures, which is closely related to the formula for Dirichlet en- 
ergy of Pinkall and Polthier fl993l . Xu [2004] studied the con- 
vergence of a number of methods for estimating the mean cur- 
vatures (and more generally, the Laplace-Beltrami operators). It 
was concluded that the cotangent formula does not converge except 
for some special cases, as also noted in I Hildebrand t et al. 20061 
|Wardetzky 2007 ] . Langer et al. 1 2005a. 2005b | proposed a tangent- 
weighted formula for estimating mean-curvature vectors, whose 
convergence also relies on special symmetric patterns of a mesh. 
Agam and Tang 1 2005 1 proposed a framework to estimate curva- 
tures of discrete surfaces based on local directional curve sampling. 
Cohen-Steiner and Morvan [2003 1 used normal cycle to analyze 
convergence of curvature measures for densely sampled surfaces. 

Vertex-based quadratic or higher-order polynomial fittings can pro- 
duce convergent normal and curvature estimations. Meek and Wal- 
ton (2000 1 studied the convergence properties of a number of es- 
timators for normals and curvatures. It was further generalized 
to higher-degree polynomial fittings by Cazals and Pouget 1 2005 1. 
These methods are most closely related to our approach. It is 
well-known that these methods may encounter numerical difficul- 
ties at low-valence vertices or special arrangements of vertices 
ILancaster and Salkauskas 1986 1, which we address in this paper. 
Razdan and Bae ] 2005 1 proposed a scheme to estimate curvatures 
using biquadratic Bezier patches. Some methods were also pro- 
posed to improve robustness of curvature estimation under noise. 
For example, Ohtake et al. 1 2004 1 estimated curvatures by fitting 
the surface implicitly with multi-level meshes. Recently. Yang et 
al. (2006] proposed to improve robustness by adapting the neigh- 
borhood sizes. These methods in general only provide curvature 
estimations that are meaningful in some average sense but do not 
necessarily guarantee convergence of pointwise estimates. 

Some of the methods estimate the second-order differential quan- 
tities from the surface normals. Goldfeather and Interrante 1 2004 1 
proposed a cubic-order formula by fitting the positions and normals 
of the surface simultaneously to estimate curvatures and principal 
directions. Theisel et al. 1 2004 1 proposed a face-based approach for 
computing shape operators using linear interpolation of normals. 
Rusinkiewicz [2004| proposed a similar face-based curvature esti- 
mator from vertex normals. These methods rely on good normal 
estimations for reliable results. Zorin and coworkers IZori n 20051 
|Grinspun et al. 2006 1 proposed to compute a shape operator using 
mid-edge normals, which resembles and "corrects" the formula of 
IQriate and Cervera 19931 . Good results were obtained in practice 
but there was no theoretical guarantee of its order of convergence. 

3 Consistent Formulas of Curvatures for 
Continuous Surfaces 

In this section, we first consider the computation of differen- 
tial quantities of continuous surfaces (as opposed to discrete sur- 
faces). This is a classical subject of differential geometry (see e.g., 
Ido Carmo 1 9761). but the differential quantities are traditionally 
expressed in terms of the first and second fundamental forms, which 
are geometrically intrinsic but sometimes difficult to understand 
when a change of coordinate system is involved. Furthermore, as 



we will show shortly, it is not straightforward to evaluate the classi- 
cal formulas consistently (i.e., in a backward-stable fashion) in the 
presence of round-off errors due to the potential loss of symmetry 
and orthogonality. Starting from these classical formulas, we derive 
some simple formulas by a novel use of the singular value decom- 
position of the Jacobian of the height functions. The formulas are 
easier to understand and to evaluate. To the best of our knowledge, 
some of our formulas (including those for the symmetric shape op- 
erator and principal curvature tensor) have not previously appeared 
in the literature. 

3.1 Classical Formulas for Height Functions 

We first review some well-known concepts and formulas in differ- 
ential geometry and geometric modeling. Given a smooth surface 
defined in the global xyz coordinate system, it can be transformed 
into a local uvw coordinate system by translation and rotation. As- 
sume both coordinate frames are orthonormal right-hand systems. 
Let the origin of the local frame be [a;o, J/o, zo] (note that for con- 
venience we treat points as column vectors). Let ti and £2 be the 
unit vectors in the global coordinate system along the positive di- 
rections of the u and v axes, respectively. Then, rh — ti x £2 is 
the unit vector along the positive w direction. Let Q denote the 
orthogonal matrix composed of column vectors ti, £2, and rh, i.e., 
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(1) 



where '[' denotes concatenation. Any point x on the surface is 
then transformed to a point [u, v, f(u)] = Q T ix — xq), where 
u = (u,v). In general, / is not a one-to-one mapping over the 
whole surface, but if the uv plane is close to the tangent plane at a 
point x, then / would be one-to-one in a neighborhood of x. We 
refer to this function f(u) : R 2 — > R (more precisely, from a subset 
of R 2 into R) as a height function at x in the uvw coordinate frame. 
In the following, all the formulas are given in the uvw coordinate 
frame unless otherwise stated. 

If the surface is smooth, the height function / defines a smooth 



surface composed of points p(u) = [u,v, f(u)] T 

fv 



v/ 



Let 



denote the gradient of / with respect to u and 
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the Hessian of /, where f uv = f vu . The 



Jacobian of p(u) with respect to u, denoted by J, is then 
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The vectors p u and p v form a basis of the tangent space of the 
surface at p, though they may not be unit vectors and may not be 

orthogonal to each other. Let du denote 



dv 



The first funda- 



mental form of the surface is given by the quadratic form 



I(du) = du Gdu, where G 
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G is called the first fundamental matrix. Its determinant is g = 
det(G) = 1 + fu + fv- We introduce a symbol i = ^fg, which 
will reoccur many times throughout this section. Geometrically, 
I = \\p u x p v ||, so it is the "area element" and is equal to the area 



of the parallelogram spanned by p u and p v . The unit normal to the 
surface is then 



with 
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The normal in the global xyz frame is then 
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j(m-futi-f v t 2 ). (4) 



The second fundamental form in the basis {p u ,p v } is given by the 
quadratic form 

II(du) = du Bdu 

where 
B = — 
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B is called the second fundamental matrix, and it measures the 
change of surface normals. It is easy to verify that 



B = H/£ 
by plugging (O and (O into {5). 



(6) 



In differential geometry, the well-known Weingarten equations read 
[n u | n„] = — [p u | p„] W, where W is the Weingarten matrix 
(a 2 x 2 matrix a.k.a. the shape operator) at point p with basis 
{.PuiPv}- left-multiplying J T on both sides, we have B — 
GW, and therefore, 



W 



= G _1 B = i(J T J) 



(7) 



The eigenvalues ki and K2 of the Weingarten matrix W are the 
principal curvatures at p. The mean curvature is Kg = (/tj + 
K2)/2 and is equal to half of the trace of W. The Gaussian curva- 
ture at p is kg = «i K2 and is equal to the determinant of W . Note 
that k% > kg, because 4(«h — kg) = («i — K2) 2 . Let di and 
d2 denote the eigenvectors of W. Then ei = Jdi|| and 

e2 = Jda/ 1| Jda || are the principal directions. 

The principal curvatures and principal directions are often used to 
construct the 3x3 matrix 



C = /tieiei 



- K 2 e 2 e2". 



(8) 



We refer to C as the principal curvature tensor in the local coordi- 
nate system. Note that C is sometimes simply called the curvature 
tensor, which is an overloaded term (see e.g., |do Carmo 1992 1), 
but we nevertheless will use it for conciseness in later discussions 
of this paper. Note that n must be in the null space of C, i.e., 
Cn = 0. The curvature tensor is particularly useful because it 
can facilitate coordinate transformation. In particular, the curvature 
tensor in the global coordinate system is C g = QCQ T . 

3.2 Consistency of Classical Formulas 

The Weingarten equation is a standard way for defining and com- 
puting the principal curvatures and principal directions in differen- 
tial geometry |do Car mo 19761 , so it is important to understand its 
stability and consistency. We notice that the singular- value decom- 
position (SVD) of the Jacobian matrix J is given by 



J = 



c/£ 
s/£ 
|V/||/i 
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and s = sin 9 
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where 8 — arctan(/„// u ). In the special case of /„ = /„ = 0, we 
have £ — 1, and any 6 leads to a valid SVD of J. To resolve this 
singularity, we define 9 — when fu~fv= 0. 

In (|9), the column vectors of U are the left singular vectors of J, 
which form an orthonormal basis of the tangent space, the diagonal 
entries of S are the singular values, and V is an orthogonal matrix 
composed of the right singular vectors of J. The condition number 
of a matrix in 2-norm is defined as the ratio between its largest and 
the smallest singular values, so we obtain the following. 

Lemma 1 The condition number of J in 2-norm is I and the con- 
dition number of the first fundamental matrix G is £ 2 . 

Note that the condition number of G is equal to its determinant. 
This is a coincidence because one of the singular values of J and 
in turn G happens to be equal to 1. Suppose H has an absolute 
error SH of C>(e| | || 2 ) - From the standard perturbation analysis, 
we know that the relative error in W computed from (|7) can be 
as large as 0(l 2 e). The mean and Gaussian curvatures, if com- 
puted from the trace and determinant of W, would then have an 
absolute error of G>(^ 2 e || 1 1 2)- Therefore, the Weingarten equa- 
tion must be avoided if I 2 3> 1, and it is desirable to keep £ as 
close to 1 as possible. Note that W is nonsymmetric, and more 
precisely it is nonnormal, i.e., 7^ WW T . The eigen- 

values of nonsymmetric matrices are not necessarily real, and the 
eigenvectors of a nonnormal matrix are not orthogonal to each other 
IGolub and Van Loan 19961 . Therefore, in the presence of round- 
off or discretization errors, the principal curvatures computed from 
the Weingarten matrix are not guaranteed to be real, and the princi- 
pal directions are not guaranteed to be orthogonal, causing potential 
inconsistencies. 

Unlike the Weingarten matrix W, a curvature tensor is always 
symmetric. The curvature tensors are important concepts and 
are widely used nowadays (see e.g., (Gatzke and Grim m~2006l 
ITheisel et al. 20041 ). In principle, the two larger eigenvalues (in 
terms of magnitude) of a curvature tensors corresponding to the 
principal curvatures and their corresponding eigenvectors give a set 
of orthogonal principal directions. As long as the estimated curva- 
tures are real, one can use l[8j to construct a curvature tensor C, 
except that C would be polluted by the errors associated with the 
nonorthogonal eigenvectors of W. This pollution can in turn lead 
to large errors in both the eigenvalues and eigenvalues of C. For 
example, for a spherical patch z = y4— (x — 0.5) 2 — (y — 0.5) 2 
over the interval [0, 1] x [0, 1], we observed a relative error of up to 
8.5% in the eigenvalues of the curvature tensors constructed from 
a nonsymmetric Weingarten matrix. Another serious problem as- 
sociated with curvature tensors is that if the minimum curvature is 
close to zero, then the eigenvectors of a curvature tensor are ex- 
tremely ill-conditioned. Therefore, one must avoid computing the 
principal directions as the eigenvectors of the curvature tensor if the 
minimum curvature is close to 0. 



3.3 Simple, Consistent Formulas 

To overcome the potential inconsistencies of the classical formu- 
las of principal curvatures and principal directions in the presence 
of large I or round-off errors, we derive a different set of explicit 
formulas. Our main idea is to enforce orthogonality and symme- 
try explicitly and to make the error propagation independent of I as 
much as possible, starting from the following simple observation. 

Lemma 2 Let U = [ui j 1x2] be a 3x2 matrix whose column vec- 
tors form an orthonormal basis of the Jacobian, and C be the cur- 



Mature tensor. U T CU is the shape operator in the basis {iti, U2}, 
and it is a symmetric matrix with an orthogonal diagonalization 
XAX~ X , whose eigenvalues (i.e., the diagonal entries of A.) are 
real and whose eigenvectors are orthonormal (i.e., X T = X -1 ). 
The column vectors ofUX are the principal directions. 

PROOF. Let Ki and di be the eigenvalues and eigenvectors of the 
Weingarten matrix W in {7}. Therefore, 



Bdi = GWdi = KiGdi. 



(11) 



Let S = U T J. By left-multiplying S T on both sides of i I 1 [ 1, we 



have 



S~ T BS~ 



l (sd t ) = ^(s^Gd,) = Ki(Sdi 



Therefore, the eigenvalues of S~ T BS 1 are the principal curva- 
tures m and its eigenvectors are the principal directions d\ = Sdi 
and d 2 = Stfe in the basis {ui, u 2 }, so S~ T BS~ X is the shape 
operator in this basis. This shape operator is symmetric, so its 
eigenvalues are real and its eigenvectors are orthonormal. By defi- 
nition of the curvature tensor, 

C = KiUdidiU T + K 2 Ud 2 d 2 U T 
= U(Kidid 1 + n 2 d 2 d 2 )U 
= US- T BS~ 1 U T , 

so U T CU is equal to the shape operator S~ T BS 1 . ■ 

Lemma |2] holds for any orthonormal basis of the Jacobian. There 
is an infinite number of such bases. Algorithmically, such a basis 
can be obtained from either the QR factorization or SVD of J. We 
choose the latter for its simplicity and elegance, and it turns out to 
be particularly revealing. 

Theorem 1 Let J = U"SV T be the reduced SVD of the Jacobian 
of a height function f as given by ||9J- Let U\ and u 2 denote the 
column vectors ofU and let S = S V T . The shape operator in the 
orthonormal basis {ui, u 2 } is the symmet 



ric matrix 



W = S^BS' 1 = - 



where c and s are defined in ilO\ 
PROOF. Note that 

S' 1 = VE" 1 = 



cji s/l 
— s c 



H 



c/£ 
s/l 



(12) 



1 _ 


c —s 




' l/l " 




" c/e -s ' 




s c 




1 




s/e c 



Following the same argument as for Lemma[2] Eq. dl2t holds. ■ 

To the best of our knowledge, the symmetric shape operator W 
defined in dl2t has not previously appeared in the literature. Be- 
cause W is symmetric, its eigenvalues (i.e., principal curvatures) 
are guaranteed to be real, and its eigenvectors (i.e., principal di- 
rections) are guaranteed to be orthonormal. In addition, because 
W is given explicitly by J 1 2b . the relative error in H is no longer 
amplified by a factor of I 2 (or even £) in the computation of W. 
Therefore, this symmetric shape operator provides a consistent way 
for computing the principal curvatures or principal directions. Fur- 
thermore, using this result we can now compute the curvature tensor 
even without forming the shape operator. 

Theorem 2 Let the transformation from the local to global coordi- 
nate system be x — Qu + xo. The curvature tensors in the local 



and global coordinate systems are 



C = J +T B J 4 



C g Jg^ BJg^ 
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< +T HJ„. 



J +1 HJ + 
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respectively, where 



1 "t" fv fufv fu 

fufv 1 H~ fu fv 



(13) 
(14) 

(15) 



is the pseudo-inverse of J, and = J + Q T is the pseudo-inverse 
of J 3 = QJ. In addition, H = £J T CJ = U T g C g J g . 

PROOF. The principal directions are ej = Ud\ and e 2 = Ud 2 in 
the local coordinate system, where di are the eigenvectors of W . 
Therefore, 

C = Kieie-i + K 2 e 2 e 2 
= UWU T 

= (W 1 V T )B(VV 1 U T ) 
= J +T BJ+. 

By left multiplying £J T and right multiplying J of both sides, we 
then have 

£J T CJ = £J T J +T BJ + J = £B = H. 

The curvature tensor in the global coordinate system is C g = 
QCQ T , which results in Qj) and H = £J T g C g J g . To obtain 
the explicit formula for J + , simply plug in c, s, and £ into 



Alternatively, we may obtain J + by expanding and simplifying 

(j t j)- 1 j 4 '. " m 

Theorem[2]allows us to transform between the curvature tensor and 
the Hessian of a height function, both of which are symmetric ma- 
trices. Furthermore, we can easily compute the gradient and Hes- 
sian of one height function from their corresponding values of an- 
other height function in a different coordinate system corresponding 
to the same point of a smooth surface. 
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Corollary 1 Given the gradient V/ = 



A. 

fv 



and Hessian H 



of a height function f in an orthonormal coordinate system, let £ 
denote %/l + /u + /». n the unit normal [— fu,— fv,l] /£, C 
the curvature tensor J +T H J + / £, where J = [I 2x2 | V/] T . Let 
Q = [q 1 I q 2 I q 3 ] be the orthogonal matrix whose column vectors 
form the axes of another coordinate system, where q :i n > 0. Let 

[a, j3, 7] T denote Q n, and let f denote the height function in the 
latter coordinate frame corresponding to the same smooth surface. 

-a/7 



Then the gradient of f is V/ 



-/J/7 



and the Hessian of f 



is J T CJ/7, where J = Q [/ 2X 2 V/ 

A direct implication of this corollary is that knowing the consis- 
tent surface normal and curvature tensor of a smooth surface in 
the global coordinate frame is equivalent to knowing the consis- 
tent gradient and Hessian of its corresponding height function in 
any local coordinate frame in which the Jacobian is nondegenerate 



(i.e., I > 0). Here, the consistency refers to the fact that the normal 
n is in the null space of the curvature tensor C and that both C 
and the Hessian H are symmetric (i.e., Cn = 0, C T = C, and 
H = H T ). 

The preceding formulas for the symmetric shape operator and 
curvature tensor appear to be new and are particularly useful 
when computing the principal curvatures and principal direc- 
tions. Because many applications require the mean curvature and 
Gaussian curvature, we also give the following simple formulas, 
which are equivalent to those in classical differential geometry 
Ido Carmo 1 976 , p. 163] but are given here in a more concise form. 

Theorem 3 The mean and Gaussian cur\>ature of the height func- 
tion /(it) : R 2 Rare 



K H = 



tr(H) 
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(V/) T /f(V/) det(H) 
■ , ana kg = 



2 p 



/ 4 



(16) 



PROOF. Let v = [c,s] T and v x = 
defined in dipt . The trace of W is 
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where the last step uses v T Hv + u XT i3"i; x = tr(JT) 
V/, therefore, 



Since 
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With regard to kg, we have 



kg = det(W) = det(5 _1 ) 2 det(i//^) 



det(iJ) 

£ 4 ' 



The same argument may be applied to other methods that assume 
more than five independent parameters per point for the first- and 
second-order differential quantities. For example, if the differen- 
tial quantities are evaluated from a parameterization of a surface 
where x, y, and z coordinates are viewed as independent functions 
of u and v, one must compute more than five parameters, so their 
numerical solutions may not be consistent. The gradient and the 
Hessian of the height function contain exactly five degrees of free- 
dom when the symmetry of the Hessian is enforced, so we have 
reduced the problem consistently into the computation of the gra- 
dient and Hessian of the height function. We therefore build our 
method based on the estimation of the gradient and Hessian of the 
height functions. 

4 Computing Gradient and Hessian of Height 
Function 

To apply the formulas for continuous surfaces to discrete 
surfaces, we must select a region of interest, define a lo- 
cal uvw coordinate frame, and approximate the gradient and 
Hessian of the resulting height function. We build our 
method based on a local polynomial fitting. Polynomial fit- 
ting is not a new idea and has been studied intensively 
(e.g., |Lancaster and Salkauskas 1986, de Boor and Ron 1992, 
IMeek and W alton 2000, Cazals and Pouget 2005]). However, it is 
well-known that polynomial fitting may suffer from numerical in- 
stabilities, which in turn can undermine convergence and lead to 
large errors. We propose some techniques to overcome instabili- 
ties and to improve the accuracy of fittings. These techniques are 
based on classical concepts in numerical linear algebra but are cus- 
tomized here for derivative computations. In addition, we propose a 
new iterative-fitting scheme and also present a unified error analysis 
of our approach. 

4.1 Local Polynomial Fitting 

The local polynomial fitting can be derived from the Taylor series 
expansion. Let f(u) denote abivariate function, where u = (u, v), 



and let Cjk be a shorthand for 



Qj + k 



pansion of / about the origin uq = (0, 0) is 



/(0). The Taylor series ex- 



This completes our description of the explicit formulas for the 
differential quantities. Since their derivations, especially for the 
principal curvatures, principal directions, and curvature tensor, are 
based on the shape operator associated with the left singular- vectors 
of the Jacobian, many intermediate terms cancel out due to symme- 
try and orthogonality, and the final formulas are remarkably sim- 
ple. It is important to note that even if the gradient or Hessian 
contain input errors, as long as the Hessian is symmetric and the 
Jacobian is nondegenerate, our formulas are consistent in the fol- 
lowing sense: The principal curvatures are real, the principal direc- 
tions are orthonormal, and the surface normal is orthogonal to the 
column space of J, the row spaces of J + , and the column space of 
the curvature tensor. 

The consistency of our formulas is enabled by a simple fact: A con- 
sistent set of first- and second-order differential quantities have five 
degrees of freedom, which is a direct result of Corollary Q] In the 
Weingarten equations, numerically there are six degrees of freedom 
(three in G and three in B), so there is an implicit constraint (which 
should have enforced the orthogonality of the principal directions) 
not present in the equations. In general, if a system of equations has 
more degrees of freedom than the intrinsic dimension of the prob- 
lem (i.e., with implicit constraints), its numerical solution would 
likely lead to inconsistencies in the presence of round-off errors. 
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Given a positive integer d, a function f(u) can be approximated to 
(d + l)st order accuracy about the origin uq as 
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assuming / has d+1 continuous derivatives. The derivatives of the 
Taylor polynomial are the same as / at uq up to degree d. 

Given a set of points sampling a small patch of a smooth surface, the 
method of local polynomial fitting approximates the Taylor poly- 
nomial by estimating Cjk from the given points. The degree of the 
polynomial, denoted by d, is called the degree of fitting. We will 
limit ourselves to relatively low degree fittings (say d < 6), because 
high-degree fittings tend to be more oscillatory and less stable. To 
estimate Cjk at a vertex of a surface mesh, we choose the vertex to 
be the origin of a local uvw coordinate frame, where the w compo- 
nent of the coordinates in this frame would be the height function /. 



We use an approximate vertex normal (e.g., obtained by averaging 
the face normals) as the w direction, so that the condition number 
of the Jacobian of / would be close to 1. Plugging in each given 
point [m,Vi, fi] into i l 1 St . we obtain an approximate equation 



is to minimize a weighted norm (or semi-norm), i.e., 
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which has n = (d + l)(d + 2)/2 unknowns (i.e., Cjk for < 
j + k < d, j > and k > 0). As a concrete example, for cubic 
fitting the equation is 
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Let m denote the number of these given points (including the vertex 
uo itself), we then obtain an m x n rectangular linear system. If 
we want to enforce the fit to pass through the vertex itself, we can 
simply set coo = and remove the equation corresponding to Mo, 
leading to an (m — 1) X (n — 1) rectangular linear system. In 
our experiments, it seems to have little benefit to do this, so we do 
not enforce coo = in the following discussions. However, our 
framework can be easily adapted to enforce coo = if desired. 

Suppose we have solved this rectangular linear system and obtained 



the approximation of Cjjt. At uq, the gradient of / is 



cio 
coi 



and the Hessian is 



Plugging their approximations 



C20 Cn 
Cll C 2 

into the formulas in Section [3] would give us the normal and cur- 
vature approximations at uo- This approach is similar to the fit- 
ting methods in |Cazals and Pouget 2005 , Meek and Walton 2000 1. 
Note that the Taylor series expansions of f u and f v about uq are 
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where the residual is 0(||u|| ). The expansions for the second 
derivatives have similar patterns. Plugging in the approximations 
of Cjfc into these series, we can also obtain the estimations of the 
gradient, Hessian, and in turn the curvatures at a point u near uq. 
The remaining questions for this approach are the theoretical issue 
of solving the rectangular system as well as the practical issues of 
the selection of points and robust implementation. 

4.2 Weighted Least Squares Formulation 

Let us denote the rectangular linear system obtained from dl9t for 
i — 1 , . . . , m as 

Vc « /, (23) 

where c is an n- vector composed of Cjh, f is an m- vector com- 
posed of fi, and V is a generalized Vandermonde matrix. For 
now, let us assume that m > n and V has full rank (i.e., its col- 
umn vectors are linearly independent). Such a system is said to be 
over-determined. We will address the robust numerical solutions 
for more general cases in the next subsection. 

The simplest solution to J23t is to minimize the 2-norm of the resid- 
ual Vc — f, i.e., mine \\Vc — f\\2- A more general formulation 
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where W is an m x m diagonal matrix with nonnegative entries. 
We refer to W as a weighting matrix. Such a formulation is called 
weighted least squares |Golub and Van Loan 1996. p. 265], and it 
has a unique solution if WV has full rank. It is equivalent to a 
linear least squares problem 



Ac » b, where A = WV and b = Wf. 



(25) 



Let u>i denote the ith diagonal entry of W. Algebraically, u>i as- 
signs a weight to each row of the linear system. Geometrically, it 
assigns a priority to each point (the larger Ui, the higher the prior- 
ity). If / is in the column space of V, then a nonsingular weighting 
matrix has no effect on the solution of the linear system. Other- 
wise, different weighting matrices may lead to different solutions. 
Furthermore, by setting u>i to be zero or close to zero, the weighting 
matrix can be used as a mechanism for filtering out outliers in the 
given points. 

The weighted least squares formulation is a general technique, but 
the choice of W is problem dependent. For local polynomial fit- 
ting, it is natural to assign lower priorities to points that are farther 
away from the origin or whose normals differ substantially from the 
w direction of the local coordinate frame. Furthermore, the choice 
of W may affect the residual and also the condition number of A. 
Let rhi denote an initial approximation of the unit normal at the 
ith vertex (e.g., obtained by averaging face normals). Based on the 
above considerations, we choose the weight of the ith vertex as 
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where 7^ = max(0, ihj mo) and e 
general, this weight is approximately equal to ||iti| 
the exponent d/2 tries to balance between accuracy and stability. 
The factor approaches 1 for fine meshes at smooth areas, but 
it serves as a safeguard against drastically changing normals for 
coarse meshes or nonsmooth areas. The term e prevents the weights 
from becoming too large at points that are too close to uq . 

The weighting matrix scales the rows of V. However, if Ui or Vi 
are close to zero, the columns of V can be poorly scaled, so that the 
ith row of A would be close to [wj, 0, . . . , 0]. Such a matrix would 
have a very large condition number. A general approach to alleviate 
this problem is to introduce a column scaling matrix S. The least 
squares problem then becomes 
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(27) 



Here, S is a nonsingular n x n diagonal matrix. Unlike the weight- 
ing matrix W, the scaling matrix S does not change the solution 
of c under exact arithmetic. However, it can significantly improve 
the condition number and in turn improve the accuracy in the pres- 
ence of round-off errors. In general, let aj denote the jth column 
vector of A. We choose S = diag(l / 1| o-i || 2 , • • ■ , l/||a„||2), as it 
approximately minimizes the 2-norm condition number of AS (see 
I Gol ub and Van Loan 19961 p. 265] and Ivan der Sluis 19691 ). 

4.3 Robust Implementation 

Our preceding discussions considered only over-determined sys- 
tems. However, even if m > n, the matrix V (and in turn AS) 
may not have full rank for certain arrangements of points, so the 
weighted least squares would be practically under-determined with 



an infinite number of solutions. Numerically, even if V has full 
rank, the condition number of AS may still be arbitrarily large, 
which in turn may lead to arbitrarily large errors in the estimated 
derivatives. We address these issues by proposing a systematic way 
for selecting the points and a safeguard for QR factorization for 
solving the least squares problem. 

4.3.1 Section of Points 

As pointed out in | Lancaster and Salkauskas 1986 1 , the condition 
number of polynomial fitting can depend on the arrangements of 
points. The situation may appear to be hopeless, because we in gen- 
eral have little control (if any at all) about the locations of the points. 
However, the problem can be much alleviated when the number of 
points m is substantially larger than the number of unknowns n. To 
the best of our knowledge, no precise relationship between the con- 
dition number and m/n has been established. However, it suffices 
here to have an intuitive understanding from the geometric inter- 
pretation of condition numbers. Observe that an m x n matrix M 
(with m > n) is rank deficient if its row vectors are co-planar (i.e, 
contained in a hyperplane of dimension less than n), and its con- 
dition number is very large if its row vectors are nearly co-planar 
(i.e., if they lie in the proximity of a hyperplane). Suppose the row 
vectors of M are random with a nearly uniform distribution, then 
the probability of the vectors being nearly co-planar decreases ex- 
ponentially as m increases. Therefore, having more points than 
the unknowns can be highly beneficial in decreasing the probability 
of ill-conditioning for well-scaled matrices. In addition, it is well- 
known that having more points also improves noise resistance of 
the fitting, which however is beyond the scope of this paper. 

Based on the above observation, we require m to be larger than n, 
but a too large m would compromise efficiency and may also under- 
mine accuracy. As a rule of thumb, it appears to be ideal for m to be 
between 1.5n and 2n. For triangular meshes, we define the follow- 
ing neighborhoods that meet this criterion. First, let us define the 
1-ring neighbor faces of a vertex to be the triangles incident on it. 
We define the 1-ring neighborhood of a vertex to be the vertices of 
its 1-ring neighbor faces, and define the 1.5-ring neighborhood as 
the vertices of all the faces that share an edge with a 1-ring neigh- 
bor face. This definition of 1-ring neighborhood is conventional, 
but our definition of 1.5-ring neighborhood appears to be new. Fur- 
thermore, for k > 1, we define the (k + l)-ring neighborhood of 
a vertex to be the vertices of the fc-ring neighborhood along with 
their 1-ring neighbors, and define the (k + 1.5)-ring neighborhood 
to be the vertices of the fc-ring neighborhood along with their 1.5- 
ring neighbors. Figure Q] illustrates these neighborhood definitions 
up to 2.5 rings. 

Observe that a typical ^ii-ring neighborhood has about the ideal 
number of points for the dth degree fittings at least up to degree 
six, which is evident from Table Q] Therefore, we use this as the 
general guideline for selecting the points. Note that occasionally 
(such as for the points near boundary) the ^±I-ring neighborhood 
may not have enough points. If the number of points is less than 
1.5n, we increase the ring level by 0.5 up to 3.5. Note that we do 
not attempt to filter out the points that are far away from the origin 
at this stage, because such a filtering is done more systematically 
through the weighting matrix in d24t . 

4.3.2 QR Factorization with Safeguard 

Our discussions about the point selection was based on a prob- 
ability argument. Therefore, ill-conditioning may still occur es- 
pecially near the boundary of a surface mesh. The standard 
approaches for addressing ill-conditioned rectangular linear sys- 
tems include SVD or QR factorization with column pivoting (see 




2.5 ring 2 ring 



Figure 1: Schematics of 1-ring, 1.5-ring, 2-ring, and 2.5-ring 
neighbor vertices of center vertex (black) in each illustration. 



Table 1: Numbers of coefficients in dth degree fittings versus num- 
bers of points in typical =±! rings. 



degree (d) 
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#coeffs. 
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10 


15 


21 


28 


#points in ring 


7 


13 


19 


31 


37 


55 



IGolub and Van Loan 19961 p. 270]. When applied to <27>, the for- 
mer approach would seek a solution that minimizes the 2-norm of 
the solution vector ||d||2 among all the feasible solutions, and the 
latter provides an efficient but less robust approximation to the for- 
mer. Neither of these approaches seems appropriate in this context, 
because they do not give higher priorities to the lower derivatives, 
which are the solutions of interest. 

Instead, we propose a safeguard for QR factorization for the local 
fittings. Let the columns of V (and in turn of A) be sorted in in- 
creasing order of the derivatives (i.e., in increasing order of j + k). 
Let the reduced QR factorization of AS be 

AS = QR, 

where Q is m x n with orthonormal column vectors and R is 
an n x n upper-triangular matrix. The 2-norm condition num- 
ber of AS is the same as that of R. To determine whether AS 
is nearly rank deficient, we estimate the condition number of R 
in 1-norm (instead of 2-norm, for better efficiency), which can 
be done efficiently for triangular matrices and is readily avail- 
able in linear algebra libraries (such as DGECON in LAPACK 
| Anderson et al. 1999]). If the condition number of R is too large 
(e.g., > 10 3 ), we decrease the degree of the fitting by removing the 
last few columns of AS that correspond to the highest derivatives. 
Note that QR factorization need not be recomputed after decreas- 
ing the degree of fitting, because it can be obtained by removing 
the corresponding columns in Q and removing the corresponding 
rows and columns in R. If the condition number is still large, we 
would further reduce the degree of fitting until the condition num- 
ber is small or the fitting becomes linear. Let Q and R denote the 
reduced matrices of Q and R, and the final solution of c is given 
by 

c = SR 1 Q T b, (28) 



where R denotes a back substitution step. Compared to SVD or 
QR with partial pivoting, the above procedure is more accurate for 
derivative estimation, as it gives higher priorities to lower deriva- 
tives, and at the same time it is more efficient than SVD. 

4.4 Iterative Fitting of Derivatives 

The polynomial fitting above uses only the coordinates of the given 
points. If accurate normals are known, it can be beneficial to take 
advantage of them. If the normals are not given a priori, we may 
estimate them first by using the polynomial fitting described earlier. 
We refer to this approach as iterative fitting. 

First, we convert the vertex normals into the gradients of the 
height function within the local uvw frame at a point. Let hi = 
[ai,/3i,7i] denote the unit normal at the ith point in the uvw 
coordinate system, and then the gradient of the height function is 

— OLi/li 



Let a jk = c (j+1)k and b jk = c ](k+\)- By plugging 



Ui, Vi, and f u {ui,Vi) ~ — Qi/7i into J2U . we obtain an equation 
for the coefficients as; similarly for f v and 6s. For example, for 
cubic fittings we obtain the equations 
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If we enforce dj,k+i = explicitly, we would obtain a sin- 

gle linear system for all the as and 6s with a reduced number of 
unknowns. Alternatively, aj k and bj k can be solved separately us- 
ing the same coefficient matrix as ( 123b and the same weighted least 
squares formulation, but two different right-hand side vectors, and 
then a,j t k+i and bj + i, k must be averaged. The latter approach com- 
promises the optimality of the solution without compromising the 
symmetry and the order of convergence. We choose the latter ap- 
proach for simplicity. 

After obtaining the coefficients as and 6s, we then obtain the Hes- 
sian of the height function at tip as 
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From the gradient and Hessian, the second-order differential quan- 
tities are then obtained using the formulas in Section [3] We sum- 
marize the overall iterative fitting algorithm as follows: 

1) Obtain initial estimation of vertex normals by averaging face nor- 
mals for the construction of local coordinate systems at vertices; 

2) For each vertex, determine ^i-ring neighbor vertices and up- 
grade neighborhood if necessary; 

3) For each vertex, transform its neighbor points into its local coor- 
dinate frame, solve for the coefficients using QR factorization with 
safeguard, and convert gradients into vertex normals; 

4) If iterative fitting is desired, for each vertex, transform normals 
of its neighbor vertices to the gradient of the height function in its 
local coordinate system to solve for Hessian; 

5) For each vertex, convert Hessian into symmetric shape operator 
to compute curvatures, principal directions, or curvature tensor. 

In the algorithm, some steps (such as the collection of neighbor 
points) may be merged into the loops in later steps, with a trade- 
off between memory and computation. Note that iterative fitting is 



optional, because we found it to be beneficial only for odd-degree 
fittings, as discussed in more detail in the next subsection. With- 
out iterative fitting, the Hessian of the height function should be 
obtained at step 3. Finally, note that we can also apply the idea of 
iterative fitting to convert the curvature tensor to the Hessian of the 
height function and then construct a fitting of the Hessian to obtain 
higher derivatives. However, since the focus of this paper is on first- 
and second-order differential quantities, we do not pursue this idea 
further. 

4.5 Error Analysis 

To complete the discussion of our framework, we must address the 
fundamental question of whether the computed differential quanti- 
ties converge as the mesh gets refined, and if so, what is the con- 
vergence rate. While it has been previously shown that polynomial 
fittings produce accurate normal and curvature estimations in noise- 
free contexts (Cazals and Pouget 2005] IMeek and Walton 20001 . 
our framework is more general and uses different formulas for com- 
puting curvatures, so it requires a more general analysis. Let h de- 
note the average edge length of the mesh. We consider the errors in 
terms of h. Our theoretical results have two parts, as summarized 
by the following two theorems. 

Theorem 4 Given a set of points in uvw coordinate frame that 
interpolate a smooth height function f or approximate f with an 
error of 0(h d+1 ) along the w direction. Assume the point distri- 
bution and the weighing matrix are independent of the mesh reso- 
lution, and the scaled matrix AS in \27\ has a bounded condition 
number. The degree-d weighted least squares fitting approximates 
c ]k toO(h d - ] - k+1 ). 

PROOF. Let c denote the vector composed of Cj k , the exact partial 
derivatives of /. Let b denote Ac. Let r = b — b, which we 
consider as a perturbation in the right-hand side of 

Ac Rib + r. 

Because the Taylor polynomial approximates / to 0(h d+1 ), and / 
is approximated to 0(h d+1 ) by the given points, each component 
of / - Vc in <23) is 0(h d+1 ). Since r = W(f - Vc) and the 
entries in W are 0(1), each component of r is 0(h 

The error in Cj k and in r are connected by the scaled matrix AS. 
Since the point distribution are independent of the mesh resolution, 
Ui — &{h) and w, = O(h). The entries in the column of A cor- 
responding to Cj k are then Q(h-' +k ), so are the 2-norm of the col- 
umn. After column scaling, the entries in AS are then 0(1), so are 
the entries in its pseudo-inverse (AS) + . The error of d in J27l i is 
then 8d = (AS) + r. Because each component of r is 0(h d+1 ), 
each component of 5d is 0( Kh d+1 ), where k is the condition num- 
ber of AS and is bounded by a constant by assumption. The er- 
ror in c is then SSd, where the scaling factor associated with Cj k 
is 0(l/h-' +k ). Therefore, the coefficient Cj k is approximated to 

o(h d - j - h+1 ). m 

Note that our weights given in l !26t appear to depend on the mesh 
resolution. However, we can rescale it by multiplying h d ^ 2 to make 
it 0(1) without changing the solution, so Theorem[4]applies to our 
weighting scheme. Under the assumptions of Theorem [4] degree- 
d fitting approximates the gradient to 0(h d ) and the Hessian to 
0(/i d_1 ) at the origin of the local frame, respectively. Using the 
gradient and Hessian estimated from our polynomial fitting, the es- 
timated differential quantities have the following property. 

Theorem 5 Given the position, gradient, and Hessian of a height 
function that are approximated to 0(h d+1 ), 0(h d ) and 0(/i d_1 ), 



respectively, a) The angle between the computed and exact normals 
is 0(h d ); b) the components of the shape operator and curvature 
tensor are approximated 0(ft d_1 ) by | U2I ) and ||73I >,- c) the Gaus- 
sian and mean curvatures are approximated to 0(/i d_1 ) by U6h 

PROOF, a) Let /„ and /„ denote the estimated derivatives of /, 
which approximate the true derivatives f u and /„ to 0(h d ). Let 
land^denote ||[— /„, — fv, 1]|| and ||[- /„,-/„, 1]||, respectively. 
Let h denote the computed unit normal [—fu, —fv, 1] T /^ and n 
the exact unit normal [— f u , —f v , 1] T /^- Therefore, 

^ A _ n-fn,-fvAf -i[-U,-fv,lf 




Figure 2: Sample unstructured meshes of sphere and torus. 



where the numerator is 0(h d ) and the denominator is 0(1). Let 9 
denote arccos h T h. Because ||n 
O(0 4 ), it then follows that 6 is 0(h d ) 



nil! = 2 



b) In d 1 2t . if fu = fv= 0, then W = H, which is approximated 
to 0(/i d_1 ). Otherwise, H is approximated to 0(h d ^ 1 ) while £, 
c, and s are approximated to 0{h d ), so the error in W is 0(/i d_1 ). 
Similarly, in equations J 1 3b and d!4t . the dominating error term is 

the 0(h d ~ 1 ) error in H, while J + and J + are approximated to 
0(h d ) by their explicit formulas, so C and C 9 are approximated 
toO^" 1 ). 

c) In J 1 6b . £ and V/ are approximated to 0(h d ) and H is approxi- 
mated to 0(fo d_1 ). Therefore, kg and kh are both approximated 
toOO^ 1 ). ■ 

The above analysis did not consider iterative fitting. Following the 
same argument, if the vertex positions and normals are both ap- 
proximated to 0(h d+1 ') and the scaled coefficient matrix has full 
rank, then the coefficients and bjk are approximated to order 
0(/i d_ -' _fe+1 ') by our iterative fitting. The error in the Hessian 
would then be 0(h d ), so are the estimated curvatures. Therefore, 
iterative fitting is potentially advantageous, given accurate normals. 
Note that none of our analyses requires any symmetry of the in- 
put data points to achieve convergence. For even-degree fittings, 
the leading term in the remainder of the Taylor series is odd de- 
gree. If the input points are perfectly symmetric, then the resid- 
ual would also exhibit some degree of symmetry, and the leading- 
order error term may cancel out as in the centered-finite-difference 
scheme. Therefore, superconvergence may be expected for even- 
degree polynomial fittings, and iterative fitting may not be able to 
further improve their accuracies. Regarding to the principal direc- 
tions, they are inherently unstable at the points where the maximum 
and minimum curvatures have similar magnitude. However, if the 
magnitudes of the principal curvatures are well separated, then the 
principal directions would also have similar convergence rates as 
the curvatures. 



5 Experimental Results 

In this section, we present some experimental results of our frame- 
work. We focus on the demonstration of accuracy and stabil- 
ity as well as the advantages of iterative fitting, the weighting 
scheme, and the safeguarded numerical solver. We do not attempt 
a thorough comparison with other methods; readers are referred to 
| Gatzke and Grimm 2006 1 for such a comparison of earlier meth- 
ods. We primarily compare our method against the baseline fitting 
methods in |Cazals and Pouget 2005 , Meek and Walton 2000 1, and 
assess them for both closed and open surfaces. 



5.1 Experiments with Closed Surfaces 

We first consider two simple closed surfaces: a sphere with unit 
radius, and a torus with inner radius 0.7 and outer radius 1.3. 
We generated the meshes using GAMBIT, a commercial software 
from Fluent Inc. Our focuses here are the convergence rates with 
and without iterative fitting as well as the effects of the weight- 
ing scheme. For convergence test, we generated four meshes for 
each surface independently of each other by setting the desired edge 
lengths to 0.1, 0.05, 0.025, and 0.0125, respectively. Figure [2] 
shows two meshes that are coarser but have similar unstructured 
connectivities as our test meshes. 

We first assess the computations of normals using fittings of degrees 
between one and six. Figure [3] shows the errors in the computed 
normals versus the "mesh refinement level." We label the plots by 
the degrees of fittings. Let v denote the total number of vertices, and 
let hi and hi denote the exact and computed unit vertex normals at 
the ith vertex. We measure the relative L2 errors in normals as 
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We compute the convergence rates as 
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and show them at the right ends of the curves. In our tests, the 
convergence rates for normals were equal to or higher than the de- 
grees of fittings. For spheres, the convergence rates of even-degree 
fittings were about one order higher than predicted, likely due to 
nearly perfect symmetry and error cancellation. 



L 2 error of normals 



L error of normals 




Figure 3:1/2 errors in computed normals for sphere ( left) and torus 
(right). 

Second, we consider the computations of curvatures. It would be 
excessive to show all the combinations, so we only show some rep- 
resentative results. Figure [4] shows the errors in the minimum and 



maximum curvatures for the sphere. Figure [5] shows the errors in 
the mean and Gaussian curvatures for the torus. In the labels, '+' 
indicates the use of iterative fitting. Let hi and hi denote the exact 
and computed quantities at the ith vertex, we measure the relative 
errors in Li norm as 
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Let d denote the degree of a fitting. In these tests, the convergence 
rates were d — 1 or higher as predicted by theory. In addition, even- 
degree fittings converged up to one order faster due to error cancel- 
lation. The converge rates for odd-degree polynomials were about 
d — 1 but were also boosted to approximately d when iterative fit- 
ting is used. Therefore, iterative fitting is effective in improving 
odd-degree fittings. In our experiments, iterative fitting did not im- 
prove even-degree fittings. 
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Figure 4: L 2 errors in computed minimum and maximum curva- 
tures for sphere. 




Figure 6: Comparisons of curvature computations with and with- 
out weighting for sphere (left) and torus ( right). 



where (x,y) £ [0, 1] x [0, 1]. Figure|7Ja-b) shows these surfaces, 
color-coded by the mean curvatures. We use two types of meshes, 
including irregular and semi-regular meshes (see Figure |7Jc-d)). 
For convergence study, we refine the irregular meshes using the 
standard one-to-four subdivision | Gallier 2000 p. 283] and refine 
the semi-regular meshes by replicating the pattern. We computed 
the "exact" differential quantities using the formulas in Section [3] 
in the global coordinate system, but performed all other computa- 
tions in local coordinate systems. For rigorousness of the tests, we 
consider both L 2 and errors. In addition, border vertices are 
included in all the error measures, posting additional challenges to 
the tests. Note that the results for vertices far away from the bound- 
ary would be qualitatively similar to those of closed surfaces. We 
primarily consider fittings of up to degree four, since higher conver- 
gence rates may require larger neighborhoods for border vertices 
(more than 3.5-ring neighbors). To limit the length of presentation, 
we report only some representative results to cover the aforemen- 
tioned different aspects. 
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Figure 5: L 2 errors in computed mean and Gaussian curvatures 
for torus. 

The preceding computations used the weighting scheme described 
in Section l4~Tl which tries to balance conditioning and error cancel- 
lation. This weighting scheme improved the results in virtually all 
of our tests. Figure [6] shows a representative comparison with and 
without weighting (as labeled by "-nw" and "-w", respectively) for 
the maximum curvatures of the sphere and torus. 

5.2 Experiments with Open Surfaces 

We now consider open surfaces, i.e. surfaces with boundary. We 
focus on the study of stabilities and the effects of boundary and 
irregular connectivities. We use two surfaces defined by the follow- 
ing functions adopted from | Xu 2004 1 : 



z = Fi(x,y) = 
z = F 2 



1.25 + cos(5.%) 



(32) 



6 + 6(3a; - l) 2 
(x,y) =exp f~ ((x - 0.5) 2 + (y- 0.5) 2 )) , (33) 




(a) Fl 



(b) F 2 . 




(c) Irregular mesh. 



(d) Semi-regular mesh. 



Figure 7: Test surfaces (a-b) and meshes (c-d). Surfaces are color- 
coded by mean curvatures. 

We first assess the errors in the computed normals for open surfaces. 
Figure[8] shows the results for F\ over irregular meshes. We label 



all the plots and convergence rates in the same way as for closed 
surfaces. Let d denote the degree of a fitting. All these fittings de- 
livered convergence rate of d or higher in L2 errors. In L x errors, 
computed as max, ||nj — nt H2, the convergence rates were d — 0.5 
or higher, close to theoretical predictions. 
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Figure 8: L2 (left) and Lao (right) errors in computed normals for 
F\ over irregular meshes. 

Second, we assess the errors in the curvatures. Figure [9] shows the 
errors in mean curvatures for F\ over irregular meshes, and Fig- 
ure [TO] shows the errors in Gaussian curvatures for F2 over semi- 
regular meshes. Let k and k denote the exact and computed quan- 
tities. We computed the L2 error using ( 13 1 b and computed the Lao 
errors as 

max \Ri — Ki I / max{ | «; | , e} , (34) 

i 

where e = 0.01 max, \ was introduced to avoid division by too 
small numbers. The convergence rates for curvatures were approx- 
imately equal to d — 1 or higher for even-degree fittings and odd- 
degree iterative fittings. 



L error of mean curvature 
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Figure 9: L2 (left) and Lao (right) in computed mean curvatures 
for F\ over irregular meshes. 



L. error of Gaussian curvature 



L error of Gaussian curvature 




Figure 10: L2 (left) and Lao (right) in computed Gaussian curva- 
tures for F2 over semi-regular meshes. 

As noted earlier, the principal directions are inherently unstable 
when the principal curvatures are roughly equal to each other (such 
as at umbilic points). In FigureQTI we show the L2 and Lao errors 
of principal directions for F\ over semi-regular meshes, where the 
errors are measured similarly as for normals. This surface is free of 



umbilic points, and the principal directions converged at compara- 
ble rates as curvatures for iterative cubic fitting and quartic fitting. 
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Figure 11: L2 (left) and Lao (right) in computed principal direc- 
tions for F\ over semi-regular meshes. 

Finally, to demonstrate the importance and effectiveness of our con- 
ditioning procedure, Figure[T2lshows comparison of computed nor- 
mals and mean curvatures with and without conditioning (as la- 
beled by "-nc" and "-c", respectively) for fittings of degrees three, 
five, and six. Here, the conditioning refers to requiring number of 
points to be 1.5 or more times of the number of unknowns as well 
as the checking of condition numbers. Without conditioning, the 
results exhibited large errors for normals and catastrophic failures 
for curvatures, due to numerical instabilities. With conditioning, 
our framework is stable for all the tests, although it did not achieve 
the optimal convergence rate for the sixth-degree fittings due to too 
small neighborhoods for the vertices near boundary. 
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Figure 12: Comparisons of errors in computed normals (left) and 
mean curvatures (right) with and without conditioning for F\ over 
irregular meshes. 



6 Discussions 

In this paper, we presented a computational framework for comput- 
ing the first- and second-order differential quantities of a surface 
mesh. This framework is based on the computation of the first- and 
second-order derivatives of a height function, which are then trans- 
formed into differential quantities of a surface in a simple and con- 
sistent manner. We proposed an iterative fitting method to compute 
the derivatives of the height function starting from the points with or 
without the surface normals, solved by weighted least squares ap- 
proximations. We improve the numerical stability by a systematic 
point-selection strategy and QR factorization with safeguard. By 
achieving both accuracy and stability, our method delivered con- 
verging estimations of the derivatives of the height function and in 
turn the differential quantities of the surfaces. 

The main focus of this paper has been on the consistent and con- 
verging computations of differential quantities. We did not address 
the robustness issues for input surface meshes with large noise and 
singularities (such as sharp ridges and corners). We have conducted 
some preliminary comparisons with other methods, which we will 



report elsewhere. One of the major motivating application of this 
work is provably accurate and stable solutions of geometric flows 
for geometric modeling and physics-based simulations. We are cur- 
rently investigating the stability of our proposed methods for such 
problems. Another future direction is to generalize our method to 
compute higher-order differential quantities. 
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